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ABSTRACT 

Two  computer  programs  were  developed  in  advanced  BASIC 
to  perform  numerical  optimization  of  a  user  supplied  design 
problem  on  the  Hewlett  Packard  9845A  desktop  computer.   An 
executive  program,  OPCON,  provides  the  interactive  link 
between  the  computer  user  and  the  DESOP  numerical  optimiza- 
tion program.   DESOP  performs  the  numerical  optimization 
using  the  sequential  unconstrained  minimization  technique 
with  an  external  penalty  function.   The  unconstrained  sub- 
problem  is  solved  using  the  Fletcher-Reeves  method  of  con- 
jugate directions,  and  using  Golden  Section  search  and 
polynomial  interpolation  in  the  one-dimensional  search. 

A  computer  subprogram,  NISCO,  was  developed  in  advanced 
BASIC  to  model  a  nonimaging  concentrating  compound  para- 
bolic trough  solar  collector.   Thermophysical,  geophysical, 
optical  and  economic  analyses  were  used  to  compute  a  life- 
cycle  fuel  savings,  for  a  design  of  stated  thermal  capacity, 
NISCO  was  coupled  to  the  OPCON/DESOP  optimization  program 
to  find  the  design  which  maximizes  the  life-cycle  fuel 
savings. 
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NOMENCLATURE 


English  Letter  Symbols 

A  -  Solar  collector  receiver  surface  area 

A,  -  One  half  the  solar  collector  aperature  area 

B  -  Cantilevered  beam  width 

C,  -  Solar  collector  depth 

C  -  Specific  heat 

P 

CR  -  Solar  collector  concentration  ratio 

E  -  Young's  modulus  of  elasticity,  or  the 

penalty  function  exponent  used  in  DESOP 

F,  -  Lower  Golden  Section  fraction 

F~  -  Upper  Golden  Section  fraction 

H  -  Cantilevered  beam  height 

H  -  Equality  constraint  vector 

Icalc  -  A  user's  flag  in  the  DESOP  program  for 

initial  and  final  user  generated  output 

L  -  Cantilevered  beam  length,  or  the  solar 

collector  length 

m  -  Solar  collector  mass  flow  rate 

Mfr  -  Solar  collector  mass  flow  rate 

Ndv  -  Number  of  design  variables 

Obj  -  Objective  function 

Opj  -  Penalized  objective  function 

P  -  Cantilevered  beam  load 

Qa  -  Solar  collector  heat  available 


Qu 
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-  Solar  collector  heat  gain 

-  Yearly  solar  collector  heat  available 

-  Heat  flux 

-  Penalty  parameter  used  in  DESOP  or  the 
solar  collector  receiver  radius 

-  Solar  collector  receiver  radius 

-  Solar  collector  aperature  cover  temperature 

-  Solar  collector  coolant  exit  temperature 

-  Solar  collector  acceptance  half  angle 

-  Solar  collector  truncation  angle 

-  Solar  collector  distance  between  the 
reflector  and  a  point  tangent  to  the  receiver 

-  Cantilevered  beam  volume 

-  Solar  collector  reflector  coordinate 

-  Design  variable  vector 

-  Solar  collector  reflector  coordinate 


Greek  Letter  Symbols 

a  -  One-dimensional  search  move  parameter 

6  -  Cantilevered  beam  deflection 


0 


0. 
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-  Solar  collector  geometry  angle  measured 
from  the  collector  centerline  (See  Fig.  6) 

-  Solar  collector  acceptance  half  angle 

-  Solar  collector  truncation  angle 

-  Cantilevered  beam  shear  stress 

-  Cantilevered  beam  bending  stress 
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I.   INTRODUCTION 

A.   BACKGROUND 

Most  engineering  design  problems  contain  several  con- 
tinuous variables  and  as  such  have  an  infinite  number  of 
solutions.   The  purpose  of  optimization  is  to  find  the  best 
possible  solution  among  the  many  potential  solutions  for 
a  given  problem  in  terms  of  some  effectiveness  or  performance 
criteria.   There  are  several  methods  of  optimization.   The 
methods  may  be  classified  as  follows: 

Analytical  methods  which  use  the  classical  techniques 
of  differential  calculus  and  the  calculus  of  variations. 

Numerical  methods  which  use  past  information  to  generate 
better  solutions  to  the  optimization  problem  by  means  of 
iterative  procedures.   Numerical  methods  can  be  used  to 
solve  problems  that  cannot  be  solved  analytically. 

Graphical  methods  which  use  the  preparation  of  a  plot 
of  the  parameter  to  be  optimized  as  a  function  of  one  or 
more  variables.   This  method  although  simple  and  easy  to 
use  becomes  unmanageable  when  there  are  three  or  more 
design  variables. 

Experimental  methods  which  use  direct  experimentation 
of  the  actual  process,  the  results  of  one  experiment  being 
used  to  decide  on  where  to  perform  the  next  experiment. 
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Case  study  methods  which  evaluate  the  results  from  a 
number  of  representative  cases,  and  choose  the  "best" 
solution.   The  "best"  solution  is  thus  not  likely  to  be 
the  optimum  solution. 

Of  the  optimization  methods,  the  numerical  method  lends 
itself  to  computerized  solution.   As  design  is  an  inter- 
active process  between  the  designer  and  the  problem,  and 
the  desktop  computer  lends  itself  towards  dedicated  inter- 
active use,  the  development  of  a  numerical  optimization 
program  for  use  on  a  desktop  computer  in  an  interactive 
mode,  is  the  objective  of  this  thesis. 

To  date  a  great  amount  of  effort  has  been  spent  de- 
veloping reliable  and  efficient  optimization  programs  for 
mainframe  computers.   These  programs  are  fairly  large  and 
complex,  requiring  a  substantial  amount  of  core  space 
during  execution.   The  size  and  complexity  of  the  programs 
has  been  the  result  of  an  attempt  to  minimize  the  amount 
of  computer  time  required  to  perform  an  optimization  and 
thus  the  cost  to  the  user.   With  the  advent  and  availability 
of  desktop  computers,  there  has  been  a  sharp  reduction  in 
the  cost  of  computer  time  to  the  computer  user.   While  the 
desktop  computer  has  far  less  core  space  than  a  mainframe 
computer,  once  the  time  factor  is  removed  from  the  numerical 
optimization  process  it  is  possible  to  put  a  small  but 
reliable  numerical  optimization  program  on  a  desktop 
computer.   A  design  problem  concerning  the  optimal  geometry 
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for  a  nonimaging  parabolic  trough  solar  collector  was 
developed  to  demonstrate  the  engineering  application  of  the 
numerical  optimization  program  developed  for  this  thesis. 

B.   SCOPE 

The  numerical  optimization  of  a  given  function  may  be 
accomplished  using  many  varied  and  different  algorithms. 
Some  of  the  more  popular  methods  used  are:   random  search, 
linear  programming,  feasible  directions,  Golden  Section, 
Newton's  method  and  sequential  unconstrained  minimization. 
A  particular  optimization  program  will  use  one  or  more  of 
these  methods  to  efficiently  and  reliably  arrive  at  the 
best  solution  to  a  particular  problem. 

For  this  thesis  two  computer  programs  were  developed 
to  perform  numerical  optimization  on  a  desktop  computer. 
The  first  program  was  developed  as  an  executive  program  to 
control  the  optimization  process.   The  executive  program 
is  named  OPCON  which  stands  for  optimizer  CONtrol  program. 
OPCON  provides  the  interactive  link  between  the  user  and 
the  program  which  performs  the  numerical  optimization. 
OPCON  allows  the  user  to  input  data,  attach  a  specified 
analysis  subprogram  to  the  numerical  optimization  program 
and  control  execution  of  the  numerical  optimization  program. 
The  second  program  developed  was  the  numerical  optimization 
program,  DESOP.   DESOP  stands  for  DEsktop  Sequential  uncon- 
strained minimization  technique  Optimization  Program. 
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DESOP  performs  the  numerical  optimization  of  the  user 
supplied  problem  using  the  sequential  unconstrained  mini- 
mization technique  with  an  external  penalty  function.   The 
unconstrained  subproblem  is  solved  using  the  Fletcher- 
Reeves  method  of  conjugate  directions,  Golden  Section 
search,  and  polynomial  interpolation. 

The  third  computer  program,  NISCO,  was  developed  to 
model  a  nonimaging  concentrating  compound  parabolic  trough 
solar  collector  using  thermophysical ,  geophysical,  optical 
and  economic  analysis  to  compute  a  life-cycle  cost  for  a 
design  with  a  stated  energy  capacity.   NISCO  stands  for 
Nonimaging  concentrating  compound  parabolic  trough  Solar 
collector. 

C.   OBJECTIVE 

The  objective  of  this  thesis  was  to  develop  a  system 
of  interactive  programs  for  the  Hewlett-Packard  9845A 
desktop  computer  which  perform  numerical  optimization,  and 
to  demonstrate  the  capability  on  the  design  of  a  nonimaging 
concentrating  compound  parabolic  trough  solar  collector. 
Three  programs  were  developed  to  meet  the  objective:   an 
executive  program,  a  numerical  optimization  program  and 
a  solar  collector  analysis  program. 

The  purposes  of  the  executive  program  OPCON  are: 
1.   To  provide  a  primary  point  of  contact  for  the  com- 
puter user  from  which  to  effect  a  numerical  optimization 

on  any  number  of  user  prepared  analysis  subprograms. 
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2.  To  provide  a  standardized,  formatted  input  for  the 
design  variables,  side  constraints  and  optimizer  control 
parameters,  which  is  recognizable  by  all  the  programs  in 
the  numerical  optimization  package  of  programs. 

3.  To  control  the  operation  of  the  different  optimiza- 
tion and  design  analysis  programs  within  the  system  through' 
a  process  of  program  overlays  which  maximizes  the  computer 
space  available  for  the  design  analysis  program. 

4.  To  develop  a  program  which  is  portable  to  different 
computers  using  an  advanced  BASIC  language. 

The  purposes  of  the  numerical  optimization  program 
DESOP  are  to  develop  an  optimization  program  that: 

1.  Is  relialbe  in  reaching  a  design  optimum,  irrespective 
of  the  starting  point. 

2.  Will  arrive  near  the  design  optimum  using  default 
values  for  the  optimizer  control  variables. 

3.  Will  allow  the  user  to  monitor  the  optimization 
process  and  to  change  the  optimizer  control  variables  to 
more  efficiently  and/or  more  accurately  reach  the  design 
optimum. 

4.  Is  portable  to  different  computers  using  an  advanced 
BASIC  language. 

The  purposes  of  the  solar  collector  program  NISCO  are 
to: 
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1.  Model  a  nonimaging  concentrating  compound  parabolic 
trough  solar  collector  using  a  system  of  thermophysical/ 
geophysical,  optical  and  economic  equations. 

2.  Arrive  at  an  optimum  design  for  a  solar  collector 
given  a  stated  average  daily  heat  gain  and  a  life-cycle 
period. 
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II.   NUMERICAL  OPTIMIZATION 

A.   THE  CONCEPT  OF  NUMERICAL  OPTIMIZATION 

Consider  the  problem  of  designing  the  cantilevered  beam 
shown  in  Figure  1.   The  design  task  may  be  broken  down  into 
three  major  parts.   First,  the  objective  of  the  design  must 
be  determined,  which  in  this  case  is  to  minimize  the  weight 
of  the  beam  required  to  support  the  concentrated  tip  load  P. 
Second,  any  physical  constraints  that  may  effect  the  design 
must  be  determined.   Thirdly,  any  limits  which  exist  on  the 
design  variables  must  be  determined.   The  design  problem 
may  then  be  reduced  to  a  system  of  equations  as  follows: 
Minimize  the  volume  (V) 
V  =  B-H-L 
.  Subject  to: 

Bending  stress  (a,  ) 

ab  =  ^ITh^  "  20000  psi 
Shear  stress  (v) 


v  ■  2T5?H  "  1000°  pSi 


Deflection  under  load  (5) 
6  =  EOT"  *  l   inch 
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P=10000 


B— ► 


Figure   1.      Cantilevered   Beam  Design   Problem 
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With  geometric  constraints  such  that: 
0.5  $   B  <  0.5 
1.0  <  H  <  20.0 
H/B  <  10.0 
At  this  time  the  following  definitions  are  introduced: 

Objective  Function  -  The  parameter  which  is  to  be 
minimized  or  maximized  during  optimization.   The  objective 
function  always  occurs  on  the  left  side  of  the  equation 
unless  it  is  also  used  as  a  design  variable.   An  objective 
function  may  be  either  linear  or  nonlinear,  implicit  or 
explicit,  but  must  be  a  function  of  the  design  variables. 

Design  Variable  -  Any  parameter  which  the  optimiza- 
tion process  is  allowed  to  change  in  order  to  improve  the 
design.   Design  variables  appear  only  on  the  right  hand 
side  of  equations  in  the  analysis  program. 

Inequality  Constraint  -  Any  parameter  which  must 
not  exceed  specified  bounds  for  the  design  to  be  acceptable. 
Constraint  functions  always  appear  on  the  left  side  of 
equations.   A  constraint  may  be  linear,  nonlinear,  implicit 
or  explicit,  but  must  be  a  function  of  one  or  more  design 
variables. 

Equality  Constraint  -  Any  parameter  which  must  equal 
a  specified  value  for  the  design  to  be  acceptable.   The 
same  rules  apply  to  equality  constraints  as  inequality 
constraints. 
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Side  Constraint  -  Any  upper  or  lower  bound  placed 
upon  a  design  variable.   Side  constraints  are  usually  not 
included  in  the  system  of  equations  that  comprise  an 
analysis  program.   Instead  they  are  usually  included  as 
part  of  the  data  input  to  the  optimization  program. 

Analysis  Code  -  The  system  of  equations  utilizing 
the  design  variables  which  are  used  to  calculate  the  objective 
function  and  the  constraints  of  a  particular  design  problem. 

The  general  optimization  problem  may  thus  be  stated 
mathematically  as: 

Find  the  set  of  design  variables  X.  where  i  =  1,2,... ,n 
which  will: 

Minimize  the  objective  function  (Obj) 

Obj  =  f(X) 
Subject  to: 

Inequality  constraints  (G) 

G.  (X)  <  0        j  =  1,2, . . . ,m 

Equality  constraints  (H) 

H. (X)  =  0         j  =  1,2, ... ,1 

Side  constraints 

X*  1  X.  <  X?     i  =  1,2,...  ,n 

ill 

Returning  to  the  cantilevered  beam  problem,  it  may  be 
stated  in  the  standard  format  as  follows: 

Let  X(l)  =  B,  X(2)  =  H,  and  Obj  =  Vol  =  B-H-L 
Then  minimize  Obj  =  Vol 
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G(i)   =  ,-.r.   -  i  _   : 

:-4     =  §  -  10  <  o 

- 


/.  ^ 
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penalized  objective  function  at  the  initial  design  may  be 
found  by  taking  the  first  partial  derivative  of  Opj  with 
respect  to  the  design  variables.   The  gradient  of  the 
penalized  objective  function  defines  the  direction  of 
steepest  ascent.   In  the  case  of  the  cantilevered  beam,  it 
is  desired  to  minimize  the  objective  function;  therefore, 
the  greatest  improvement  in  design  may  be  achieved  by  moving 
in  the  negative  gradient,  or  steepest  descent  direction. 
From  the  initial  design  point  a  search  is  performed  in  the 
steepest  descent  direction  for  the  minimum  value  of  Opj  in 
that  direction.   At  the  new  minimum,  the  gradient  of  the 
penalized  objective  function  is  again  determined  and  a 
search  is  performed  in  a  conjugate  direction  until  a  second 
minimum  is  found.   Successive  iterations  are  performed 
until  the  gradient  is  found  to  be  zero  or  each  successive 
iteration  produces  a  sufficiently  small  change  in  Opj  such 
that  for  all  practical  purposes  the  minimum  has  been  found. 
At  this  time  the  penalty  function  is  increased.   If  the 
design  is  in  a  region  where  there  are  no  constraints 
violated  an  increase  in  the  penalty  function  will  not 
change  the  value  of  Opj .   If  on  the  other  hand  the  design 
is  in  an  infeasible  region  where  there  are  one  or  more 
constraints  violated,  Opj  will  be  increased,  and  the 
search  for  a  new  minimum  will  commence.   If  the  minimum 
of  the  objective  function  exists  in  the  infeasible  region, 
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the  minimum  value  for  the  objective  function  in  the  feasible 
region  will  be  approached  from  the  infeasible  region  as  the 
penalty  function  is  increased.   The  design  improvement 
process  will  terminate  when  a  zero  gradient  is  found  or 
successive  iterations  produce  a  sufficiently  small  change 
in  the  value  of  Opj  and  an  increase  in  the  penalty  function 
causes  no  change  in  Opj . 

The  minimum  thus  reached  by  the  optimization  process  is 
a  minimum  with  respect  to  the  penalized  objective  surface 
immediately  surrounding  the  final  design  point.   The  opti- 
mization process  cannot  distinguish  between  local  and  global 
minimum  points.   It  is  thus  good  engineering  practice  to 
run  several  optimizations  for  a  particular  design  problem 
from  several  different  initial  design  points.   If  optimiza- 
tions performed  from  different  initial  design  points  converge 
on  the  same  minimum  point,  that  point  is  probably  a  global 
minimum.   If  on  the  other  hand  two  or  more  minimums  are 
found,  there  may  be  local  minimums  located  in  the  design 
space  being  considered  and  care  must  be  taken  to  find  the 
global  minimum. 

B.   THE  DESOP  NUMERICAL  OPTIMIZATION  PROGRAM 

The  DEsktop  Sequential  unconstrained  minimization 
technique  Optimization  Program  was  developed  using  the 
basic  optimization  approach  outlined  in  Section  II. A. 
A  copy  of  the  program  is  included  in  Appendix  E.   The  major 
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program  structure  is  shown  in  Figure  2.   The  following 
discussion  will  refer  to  Figure  2  and  describe  the  major 
features  of  the  program. 

1.   Basic  Program  Execution 

The  DESOP  program  begins  execution  when  it  is 
loaded,  linked  to  the  user's  analysis  subprogram,  and  the 
program  is  instructed  to  run  by  the  OPCON  program.   The 
above  actions  are  automatically  performed  by  the  OPCON 
program.   DESOP  is  loaded  into  the  computer  by  an  overlay 
process.   Therefore  no  variables  can  be  directly  transferred 
between  the  DESOP  and  OPCON  programs.   DESOP  begins  execu- 
tion by  reading  the  optimizer  control  variables  and  the 
design  variables  that  were  input  using  the  OPCON  program 
and  saved  to  a  mass  storage  device.   The  program  then  sets 
Icalc  equal  to  one  and  evaluates  the  objective  function 
and  constraints  at  the  initial  design  point.   Icalc  is  a 
flag  provided  the  user  to  key  user  specified  output  on  the 
initial  and  final  design  analysis.   DESOP  will  provide  the 
user  with  a  hard  copy  output  of  the  initial  design  variables, 
the  value  of  the  constraints,  the  objective  function  and 
the  penalized  objective  function.   The  user  then  has  the 
option  of  continuing  with  the  DESOP  program  to  optimize 
his  analysis  subprogram  or  to  return  to  the  OPCON  program 
to  change  one  or  more  of  the  input  parameters. 


25 


INPUT    X0 


EVALUATE    Obj     &    G 


I 


FIND  THE  GRADIENTS 
OF  Obj  AT  X 


CALCULATE  THE 
STEEPEST  DESCENT 
SEARCH  DIRECTION 


CALCULATE  THE 
FLETCHER  REEVES 
SEARCH  DIRECTION 


ESTIMATE  a  FOR  MIN.  Obj 
IN  1-D  SEARCH 


FIND  a  FOR  MIN.  Obj 
BY  EXPANDING  SEARCH, 
METHOD  OF  GOLDEN  SECTION, 
&  POLYNOMIAL  INTERPOLATION 


MOVE  TO  NEW  X 


CHECK  CONVERGENCE  OF  A Obj 


no 


INCREASE  PENALTY  FUNCTION 
RE  CHECK  CONVERGENCE 


IF  MET 


no 


END 


Figure  2.   DESOP  Flow  Diagram 
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Proceeding  with  the  optimization,  there  are  two 
major  loops  in  the  optimization  program.   The  outer  loop 
increases  the  penalty  function  when  the  inner  loop's  con- 
vergence criteria  have  been  met.   A  convergence  test  is 
then  performed  by  the  outer  loop.   If  the  convergence 
criteria  is  met,  the  optimization  process  is  considered 
finished.   If  the  convergence  criteria  for  the  outer  loop 
is  not  met,  program  execution  is  returned  to  the  inner 
loop.   The  inner  loop  performs  successive  iterations 
searching  for  the  minimum  of  the  penalized  objective 
function  with  no  increase  in  the  penalty  function  taking 
place.   When  the  inner  loop's  convergence  criteria  have 
been  met  program  execution  is  transferred  to  the  outer  loop. 

Execution  of  the  program  while  in  the  inner  loop 
proceeds  as  follows:   First,  the  gradient  of  the  penalized 
objective  function  is  calculated  by  subroutine  GRAD.   The 
program  then  computes  a  search  direction  using  either  a 
steepest  descent  method  or  the  method  of  conjugate  directions 
developed  by  Fletcher  and  Reeves  [Ref.  1].   Once  a  search 
direction  is  established  the  optimizer  attempts  to  locate 
the  minimum  value  of  the  penalized  objective  function  in 
the  search  direction.   This  process  if  referred  to  as  the 
one-dimensional  search  and  is  illustrated  in  Figure  3. 
The  efficiency  and  accuracy  to  which  the  one-dimensional 
search  for  the  minimum  of  the  penalized  objective  function 


27 


w 
o 
<D 
O 

0 

u 

Cm 

A 
O 

u 
<a 

en 


(0 

c 
o 

•H 

CO 

e 
<u 

e 

•H 

Q 
I 

a 

c 
o 

<u 

J3 
En 


m 

a 

M 

3 

•H 

fc. 


28 


is  accomplished,  is  the  key  to  successful  sequential  un- 
constrained minimization  technique  numerical  optimization. 
The  one-dimensional  problem  may  be  expressed  in  terms  of 
the  penalized  objective  function,  Opj ,  and  the  amount  of 
movement,  a,  in  the  search  direction,  S.   First  the  slope 
of  Opj  with  respect  to  a  is  calculated.   An  initial  "guess" 
of  how  far  to  move  is  made  using  subroutine  ALPGES.   The  a 
which  corresponds  to  the  minimum  value  of  Opj  in  the  one- 
dimensional  search  is  then  calculated  using  subroutine 
ALPBND  and  subroutine  QUEBIC.   The  minimum  value  of  Opj 
thus  found  is  then  compared  to  the  previous  value  of  Opj 
for  convergence  using  subroutine  CONVRG.   If  convergence 
is  not  met,  execution  returns  to  the  start  of  the  inner 
loop.   If  convergence  is  met,  execution  returns  to  the 
outer  loop. 

When  the  convergence  criteria  have  been  met  for 
both  the  inner  and  outer  loops,  the  program  proceeds  to 
set  Icalc  to  three  as  a  flag  for  user  generated  output  for 
the  final  design.   DESOP  then  provides  the  user  with  a 
hard  copy  output  of  the  final  design  variables,  objective 
function  value,  penalized  objective  value,  constraint 
values,  the  number  of  inner  loop  iterations,  the  number  of 
times  the  analysis  subprogram  was  called  and  the  final  value 
of  the  penalty  function.   The  OPCON  program  is  then  overlayed 
over  the  DESOP  program  and  program  execution  is  returned  to 
the  OPCON  program. 
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2.   Finding  the  Search  Direction 

The  first  step  in  finding  the  search  direction,  S, 
is  to  determine  the  slope  of  Opj  at  the  present  design 
point.   The  forward  finite  difference  method  is  used. where: 

a_     F(X.  +    LX.)    -   F(X. ) 

3F   _     l 3/ l   _   _ 

3X.  AX.  i 

l  l 

i  =  1,2,... Ndv 
As  3F/3X.  gives  the  direction  of  positive  slope,  the  search 
direction  is  the  negative  of  3F/3X. .   The  first  search  is 
performed  using  the  steepest  descent  as  found  above  using 
the  following  relation: 

Xl  =  X.  +  aS. 
li     i 

where  alpha  is  the  distance  moved  in  the  S  direction.   When 
a  minimum  is  obtained  along  the  direction  of  steepest 
descent,  a  new  Fletcher-Reeves  conjugate  search  direction 
[Ref.  1]  is  calculated  at  the  new  Design  point  using  the 
following  relations: 

S'  =  sf  +    BS. 

l    3X.      l 

l 


Ndv 


i 


3F  t  2 


B  =  1=1 


Ndv   r_3F_i  2 

E    L  X.  J 
i=l 


where  the  prime  denotes  values  for  the  present  iteration 
and  the  non-prime  variables  indicate  values  for  the  previous 
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iteration.   A  one-dimensional  search  is  then  performed  in 
the  new  search  direction.   Searches  are  continued  using 
the  conjugate  direction  method  for  Ndv  +  1  iterations, 
where  Ndv  is  the  number  of  design  variables.   The  search 
process  is  then  restarted  using  the  steepest  descent  method. 
The  reason  for  incorporating  the  conjugate  search  method 
is  that  the  steepest  descent  method  when  traversing  a  design 
surface  with  a  curved  valley  will  tend  to  zigzag  from  one 
side  of  the  design  surface  valley  to  the  other  making  very 
little  progress  as  is  illustrated  in  Figure  4.   The  conju- 
gate direction  method  is  much  more  efficient  in  traversing 
such  a  design  surface.   However,  as  the  conjugate  direction 
method  is  additive  upon  previous  searches,  it  has  a  tendency 
to  decrease  in  effectiveness  with  each  successive  search 
owing  to  the  accumulation  of  numerical  "noise."   For  that 
reason  the  search  process  is  restarted  with  the  steepest 
descent  method  every  Ndv  +  1  iterations,  or  when  the  con- 
jugate direction  predicts  a  positive  slope.   The  search 
direction  is  normalized  to  avoid  inaccuracies  caused  by 
numerical  ill-conditioning. 

3.   Estimating  an  Initial  Value  for  Alpha 

The  initial  estimate  for  alpha  is  made  in  the 
following  manner:   First,  the  slope  of  Opj  in  the  search 
direction  is  calculated  as  the  sum  of  each  of  the  products 
of  the  gradients  times  the  search  direction.   Then  the 
slope  of  Opj  in  the  search  direction  is  divided  by  the  value 
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of  Opj .   This  value  is  then  multiplied  by  an  improvement 
percentage  in  Opj .   This  first  estimate  is  then  applied  to 
a  series  of  conditional  tests  to  determine  the  validity  of 
the  estimate  with  respect  to  the  slope  of  Opj  and  the  magni- 
tude of  the  design  variables.   Lastly,  the  estimate  for 
alpha  is  checked  to  see  if  it  violates  any  side  constraints. 
If  it  does,  the  value  of  the  estimate  for  alpha  is  reduced 
until  the  side  constraints  are  no  longer  violated. 

4.   Calculating  Alpha 

The  calculation  of  alpha  is  the  most  critical  al- 
gorithm in  the  DESOP  program  in  providing  reliable  optimizer 
operation.   The  ability  to  accurately  and  efficiently  find 
the  minimum  of  the  penalized  objective  function  in  the  one- 
dimensional  search  affects  directly  the  operation  of  the 
optimizer.   Figure  5  illustrates  the  zigzag  phenomenon 
which  occurs  if  alpha  is  not  accurately  found.   The  zigzag 
phenomenon  is  caused  by  the  fact  that  the  optimizer  in  per- 
forming the  forward  finite  difference  for  calculating  the 
search  direction  perturbs  the  design  vector  a  very  small 
amount.   As  such  the  optimizer  can  only  "see"  the  design 
surface  that  is  immediately  adjacent  to  the  design  point. 
Therefore,  if  the  minimum  is  not  found  in  the  one-dimensional 
search,  the  optimizer  will  converge  very  slowly  on  the 
minimum. 

There  are  two  major  sections  to  the  ALPBND  subroutine. 
The  first  section  attempts  to  find  the  minimum  value  of  Opj 
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Figure  5.   The  Zigzag  Phenomenon 
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using  an  expanding  search  technique.   The  first  move  is 
the  amount  predicted  by  the  ALPGES  subroutine.   If  the 
minimum  is  not  bracketed  by  the  first  move,  additional 
moves  are  made.   Each  move  is  larger  than  the  last.   The 
size  of  the  move  is  increased  each  time  by  an  amount  equal 
to  the  size  of  the  last  move  divided  by  the  lower  Golden 
Section  fraction,  where  the  Golden  Section  fractions  are: 

Fl 2 

F  =  d  - 1 

r2      2 

The  lower  Golden  Section  fraction,  F. ,  is  used  so  that  the 
interval  will  be  consistent  with  the  Golden  Section  search 
in  the  second  section  of  the  ALPBND  subroutine.   The  expand- 
ing search  is  continued  until  the  minimum  value  of  the  ob- 
jective function  has  been  bracketed. 

Once  the  minimum  is  bracketed,  a  Golden  Section 
search  is  performed  to  reduce  the  bracketing  interval  on 
the  minimum  by  an  amount  such  that  when  the  two  end  points 
of  the  interval  are  taken  with  two  points  internal  to  the 
interval  and  a  cubic  is  passed  through  the  four  points,  the 
cubic  will  accurately  predict  the  minimum  of  the  penalized 
objective  function.   Himmelblau  in  [Ref.  2]  states  that  the 
Golden  Section  search  method  of  reducing  the  interval 
around  the  minimum  of  Opj  is  the  most  effective  of  the 
reducing  techniques  studied.   Golden  Section  search  is  based 
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on  the  splitting  of  a  line  into  two  segments  known  in 
ancient  times  as  the  "Golden  Section. "   The  ratio  of  the 
whole  line  to  the  larger  segment  is  the  same  as  the  ratio 
of  the  larger  segment  to  the  smaller  segment.   The  two 
Golden  Section  fractions  are  employed  to  split  the  interval 
bracketing  the  minimum  as  shown  in  Figure  3.   Once  the 
interval  has  been  split,  the  two  values  of  Opj  corresponding 
to  the  internal  points  are  compared  to  find  the  larger  of 
the  two.   The  internal  point  with  the  larger  value  of  Opj 
will  become  the  new  end  point  for  the  interval,  the  remain- 
ing interior  point  will  by  the  fact  that  it  was  determined 
by  a  Golden  Section  fraction,  be  equal  to  the  point  deter- 
mined by  the  other  Golden  Section  fraction.   Thus,  only  one 
new  point  must  be  calculated  to  continue  the  Golden  Section 
search.   The  search  is  continued  in  this  manner  until  the 
vertical  separation  of  the  two  end  points  with  respect  to 
the  interior  points  is  less  than  one  percent.   The  four 
values  of  the  penalized  objective  function  corresponding  to 
the  four  Golden  Section  search  points  are  then  sent  to  a 
cubic  interpolator.   The  cubic  interpolator  will  return  a 
value  for  alpha  to  predict  the  minimum  of  the  penalized 
objective  function,  and  the  minimum  of  the  cubic  function 
that  the  interpolator  has  created.   The  subroutine  ALPBND 
will  then  test  the  predicted  minimum  with  the  minimum  found 
at  the  predicted  alpha.   If  there  is  less  than  a  tenth  of 
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one  percent  difference  between  the  two  values  of  the  objec- 
tive function,  the  point  predicted  by  the  cubic  interpolator 
will  be  accepted  as  the  minimum  and  program  execution  will 
return  to  the  main  program.   If  the  predicted  minimum  is  not 
sufficiently  close  to  the  minimum  at  the  predicted  alpha, 
another  Golden  Section  search  will  be  performed  to  reduce 
the  interval  and  better  localize  the  minimum.   The  four 
points  from  the  reduced  interval  will  then  be  sent  to  the 
cubic  interpolation  subroutine.   This  process  will  continue 
until  either  the  test  for  the  minimum  is  positive  or  the 
interval  has  been  reduced  to  less  than  1E-12.   Program  exe- 
cution will  then  return  to  the  main  program. 
5.   Subroutine  QUEBIC 

Subroutine  QUEBIC  is  used  to  estimate  the  alpha  at 
which  Opj  is  a  minimum  based  on  four  point  cubic  interpola- 
tion.  If  the  function  more  closely  resembles  a  quadratic 
than  a  cubic,  a  three  point  quadratic  interpolation  is 
performed  using  the  three  points  which  bracket  the  minimum. 
If  the  predicted  minimum  is  outside  the  interval  spanned 
by  the  two  end  points  again  a  quadratic  interpolation  is 
performed.   If  the  minimum  still  lies  outside  the  two  end 
points,  the  analysis  returns  to  subroutine  ALPBND,  the 
inverval  bracketing  the  minimum  is  reduced,  and  program 
execution  returns  to  QUEBIC. 
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6.  Convergence  of  the  Penalized  Objective  Function 
The  penalized  objective  function  is  tested  for  con- 
vergence at  the  end  of  each  inner  loop  and  again  at  the  end 
of  the  outer  loop  in  the  main  program.   Convergence  is 
tested  by  calling  subroutine  CONVRG.   There  are  two  criteria 
used  for  testing  for  convergence.   The  first  tests  the  rela- 
tive difference  of  the  value  of  Opj  from  the  present  iteration 
with  the  value  of  Opj  from  the  last  iteration.   The  second 
method  tests  the  absolute  difference  of  the  two  values.   The 
second  method  is  employed  for  cases  when  the  value  of  the 
penalized  objective  function  approaches  zero.   When  conver- 
gence has  been  met  on  two  successive  iterations,  the  penalty 
function  is  increased  by  an  amount  specified  by  the  user  in 
the  executive  OPCON  program.   The  penalized  objective  function 
is  again  tested  for  convergence.   If  convergence  is  still 
met,  the  optimizer  considers  the  present  value  of  the  pena- 
lized objective  function  to  be  a  minimum,  noting  again  that 
numerical  optimization  programs  cannot  differentiate  between 
local  and  global  minimums. 

7 .  The  Penalty  Function 

The  purpose  of  the  penalty  function  is  to  increase 
the  value  of  the  objective  function  when  the  design  is  in 
an  infeasible  region.   The  infeasible  region  is  that  region 
where  one  or  more  design  constraints  are  violated.   When  a 
constraint  is  violated,  the  value  of  the  particular  constraint, 
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G.,  is  positive.   The  objective  function  is  then  penalized 
as  follows: 

Opj  =  Obj  +  R-G.E 
where : 


R  -  a  multiplication  constant 

E  -  an  exponent  constant 
This  type  of  penalty  function,  one  where  the  penalty  is 
applied  after  the  design  leaves  the  feasible  region,  is 
known  as  an  exterior  penalty  function.   The  exterior  type 
of  penalty  function  was  chosen  over  other  types,  such  as 
the  interior  or  extended  interior  penalty  function.   If  a 
function  is  discontinuous  within  the  design  space  being 
studied,  numerical  difficulties  may  be  encountered  which 
make  performing  an  optimization  of  the  design  difficult. 
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III.   SOLAR  COLLECTOR  OPTIMIZATION 

A.   THE  CONCEPT  OF  NONIMAGING  SOLAR  COLLECTORS 

At  the  present  time  there  are  numerous  schemes  for 
the  collection  of  solar  energy  and  its  conversion  to  a  more 
useful  form  of  energy.   In  the  field  of  solar  collectors 
there  are  three  broad  categories:   flat  plate  collectors, 
focusing  collectors  and  nonimaging  collectors.   The  advan- 
tages and  disadvantages  of  each  type  are  shown  in  Table  I. 
Welford  and  Winston  in  [Ref.  3]  report  that  "in  the  mid- 
19601  s,  it  was  realized  in  at  least  three  different  labora- 
tories that  light  could  be  collected  and  concentrated  for 
many  purposes,  including  solar  energy,  more  efficiently  by 
nonimaging  optical  systems  than  by  conventional  image 
forming  systems.   The  methodology  of  designing  optimized 
nonimaging  systems  differs  radically  from  conventional 
optical  design.   The  new  collectors  approach  very  closely 
the  maximum  theoretical  concentration;  and  for  two-dimensional 
geometry,  which  is  important  for  solar  energy  collection, 
this  limit  is  actually  reached."    Figure  6  shows  the  basic 
geometry  for  the  nonimaging  compound  parabolic  concentrating 


Welford,  W.  T.  and  Winston,  R. ,  The  Optics  of  Nonimaging 
Concentrators,  Light  and  Solar  Energy,  p.  ix,  Academic  Press, 
1978. 
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Figure  6.   Nonimaging  Concentrating  Solar  Collector 
Geometry 
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The  collector  depth,  Cd,  is  used  in  the  economic  analysis 
of  collector  cost.   The  truncation  angle,  9  ,  is  the  maximum 
9  used  in  determining  the  collector  geometry  for  the  trun- 
cated collector.   The  collector  may  be  significantly  trun- 
cated before  any  appreciable  change  in  the  concentration 
ratio  is  affected.   This  allows  a  savings  in  manufacturing 
costs  with  little  degradation  in  collector  performance. 

The  nonimaging  concentrating  solar  collector  will 
accept  and  deliver  to  the  absorber  all  incident  radiation 
that  falls  on  the  collector  aperature  and  that  is  within 
the  maximum  acceptance  half  angles,  6..   That  is,  there  is 
an  arc  of  sky,  29.,  from  which  all  radiation,  direct,  diffuse 
and  reflected,  is  delivered  to  the  absorber.   It  is  this 
fact  which  makes  the  nonimaging  concentrating  collector 
attractive  for  solar  energy  use.   Figures  7,  8,  and  9  show 
ray  paths  for  a  nonimaging,  truncated  collector  for  the 
following  cases:   In  Figure  7  the  solar  altitude  is  within 
the  acceptance  half  angles.   In  Figure  8  the  solar  altitude 
is  equal  to  the  acceptance  half  angle.   In  Figure  9  the 
solar  angle  is  less  than  the  acceptance  half  angle.   In 
Figure  7,  all  the  radiation  that  enters  the  collector  is 
delivered  to  the  absorber  tube  and  is  somewhat  scattered 
over  the  absorber  surface.   In  Figure  8  again  all  the 
radiation  that  enters  the  collector  is  delivered  to  the 
absorber  tube,  but  is  now  tangent  to  the  tube  and  is  con- 
centrated on  the  front  edge  of  the  absorber  tube.   In 
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Figure  7. 


Ray  Path  Drawing  for  Solar  Altitude  Within  thP 
Maximum  Acceptance  Half  Angle       Wlt^n  the 
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Figure  8.   Ray  Path  Drawing  for  Solar  Altitude  Equal  to 
the  Maximum  Acceptance  Half  Angle 
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Figure  9.   Ray  Path  Drawing  for  Solar  Altitude  Outside  the 
Maximum  Acceptance  Half  Angle 
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Figure  9  all  the  reflected  radiation  leaves  the  collector, 
thus  the  acceptance  half  angle  provides  a  very  sharp  cutoff 
angle  for  accepting  incident  radiation. 

B.   SOLAR  COLLECTOR  DESIGN  PROGRAM 

A  Nonimaging  compound  parabolic  trough  Solar  collector 
hereafter  referred  to  as  NISCO  was  chosen  to  model  in  the 
analysis  program.   There  are  seven  major  heat  flow  paths 
considered  in  the  program.   The  major  heat  flow  paths  are 
shown  in  Figure  10  and  are  outlined  below: 

q,  -  The  sum  of  the  direct,  diffuse  and  ground  reflected 
radiation  that  is  incident  on  the  collector  cover. 

q~  -  The  sum  of  the  direct,  diffuse  and  ground  reflected 
radiation  that  is  reflected  by  the  collector  cover. 

q3  -  The  sum  of  the  direct,  diffuse  and  ground  reflected 
radiation  that  is  absorbed  by  the  collector  cover. 

q.  -  The  sum  of  the  direct,  diffuse,  and  ground  reflected 
radiation  that  is  transmitted  by  the  cover  and  delivered 
either  directly  or  indirectly  to  the  absorber  and  is  absorbed 
by  the  absorber. 

q5  -  The  sum  of  the  radiation  reflected  by  the  absorber 
that  is  absorbed  by  the  cover  and  the  thermal  radiative  and 
convective  exchanges  between  the  absorber  and  the  cover. 

qfi  -  The  sum  of  the  thermal  radiative  and  the  convective 
exchanges  between  the  cover  and  the  environment. 

q_  -  The  energy  delivered  to  the  collector  coolant. 
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II 


Figure  10.   Basic  Solar  Collector  Heat  Paths 
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To  determine  the  energy  balance  on  the  collector,  a 
procedure  outlined  by  Kreith  and  Kreider,  [Ref.  6]  is  followed. 
All  important  heat  fluxes  are  first  calculated  from  basic 
heat-transfer  principles.   The  fluxes  are  then  combined  in 
heat-balance  equations  for  the  receiver,  the  aperature  cover, 
and  the  coolant  fluid.   Since  the  various  flux  terms  are 
nonlinear  in  temperature  a  simultaneous  iterative  solution 
is  used  to  solve  the  equations.   Note  that  the  first  page  of 
Appendix  G  is  a  cross  reference  list  of  the  major  equations 
used  in  the  subprogram  NISCO  and  the  sources  used  to  obtain 
the  equations. 

The  following  discussion  will  detail  the  procedure  used 
in  the  NISCO  subprogram  to  calculate  the  heat  gain  for  a 
particular  solar  collector  design  and  the  resultant  life- 
cycle  fuel  savings.   Variable  names  and  program  line  numbers 
correspond  to  those  used  in  the  NISCO  subprogram.   The  design 
variables  chosen  for  the  NISCO  subprogram  were:   (1)   Thetai, 
(9.),  the  maximum  acceptance  half  angle,  (2)  Thetat,  (9.), 
the  truncation  angle,  (3)  R,  (r) ,  the  receiver  radius,  (4)  L, 
the  collector  length,  and  (5)  Mfr,  (m) ,  the  coolant  mass  flow 
rate.   The  objective  function  of  the  NISCO  subprogram  is  the 
life-cycle  fuel  savings  of  the  collector.   Life-cycle  fuel 
savings  are  calculated  as  the  cost  savings  realized  by  the 
collector  over  purchasing  natural  gas  per  unit  quantity  of 
106  Btu's  expressed  in  present  worth  using  the  present  worth 
analysis  described  by  Newnan  in  [Ref.  13] .   The  constraints 
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placed  upon  the  design  by  the  NISCO  subprogram  were:   the 
maximum  and  minimum  values  allowed  for  the  truncation  angle 
due  to  collector  geometry  as  specified  in  Section  III. A. , 
the  mass  flow  rate  of  the  coolant  must  be  positive,  the 
receiver  radius  must  be  positive,  a  minimum  average  daily 
heat  gain,  and  a  maximum  allowable  coolant  temperature. 
The  side  constraints  placed  upon  the  design  to  ensure  that 
the  final  design  was  a  reasonable  design  were:   the  maximum 
acceptance  half  angle  must  be  greater  than  three  degrees  but 
less  than  85  degrees,  the  truncation  angle  must  be  greater 
than  185  degrees  but  less  than  260  degrees,  the  receiver 
radius  must  be  greater  than  one  tenth  inch  but  less  than  two 
inches,  the  collector  length  must  be  greater  than  five  feet 
but  less  than  100  feet,  and  the  mass  coolant  flow  rate  must 
be  greater  than  five  lbm/hr  and  less  than  1000  lbm/hr. 
The  NISCO  analysis  subprogram  proceeds  as  follows. 
A  minimum  average  daily  heat  gain  is  specified  and  entered 
as  Ql  in  line  1655.   The  subprogram  is  then  SAVED  and  the 
OPCON/DESOP  optimization  program  run.   When  DESOP  calls 
the  NISCO  subprogram,  it  will  pass  the  design  variable 
vector  X  to  the  NISCO  subprogram.   The  NISCO  subprogram  in 
lines  205  to  225  sets  the  design  variable  vector  X  equal 
to  the  design  variables  used  in  the  NISCO  subprogram. 
NISCO  then  proceeds  to  read  in  the  constants  and  data  used 
in  the  solar  collector  design.   The  design  specifications 
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are  summarized  in  Appendix  B.   Collector  geometry  calcula- 
tions are  then  performed  in  lines  82  0  to  855  to  ascertain 
the  optical  properties  of  the  collector  as  developed  in 
[Ref .  3] .   An  initial  receiver  and  aperature  cover  tempera- 
ture is  assumed  in  lines  865  and  870.   Monthly  calculations 
are  then  performed  to  calculate  the  collector  tilt  angle, 
the  minimum  accepted  solar  altitude,  the  ground  angle  factor 
and  sky  heat  loss  constants  as  specified  in  [Refs.  5  and  6]. 
An  hourly  calculation  is  then  performed  to  calculate  the 
angle  that  the  sun  makes  with  the  collector  aperature  cover, 
the  amount  of  radiation  incident  on  the  collector  aperature 
cover  and  the  collector  cover  transmissivity  and  absorptance 
as  specified  in  [Refs.  5  and  6] .    The  iterative  portion  of 
the  analysis  then  proceeds  as  follows:   First,  the  heat 
transfer  coefficients  for  the  collector  are  calculated  as 
prescribed  in  [Refs.  4,  6  and  11],  lines  1060  to  1145.   Next 
a  heat  balance  is  performed  on  the  cover  as  prescribed  in 
[Ref.  6]  and  a  new  collector  aperature  cover,  Tap,  is  cal- 
culated in  lines  1155  to  1210.   A  heat  balance  is  then 
performed  on  the  receiver  as  specified  in  [Ref.  6]  and  the 
energy  passed  to  the  coolant,  Qu,  is  calculated  in  lines 
1225  to  1305.   Finally  a  heat  balance  is  performed  on  the 
coolant  as  prescribed  in  [Ref.  6]  and  the  coolant  exit 
temperature  is  calculated  in  lines  1320  to  1365.   Knowing 
the  coolant  exit  temperature,  Tc2 ,  an  average  receiver 
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temperature  may  be  calculated  and  compared  to  the  initial 
assumed  receiver  temperature,  lines  1375  to  1405.   If  the 
new  receiver  temperature  is  within  a  tenth  of  one  percent 
of  the  old  receiver  temperature,  the  analysis  program 
proceeds  to  calculate  the  pumping  power  required  for  the 
collector.   If  the  new  receiver  temperature  is  not  within 
the  convergence  specified,  the  iterative  heat  balance 
process  is  repeated,  setting  the  aperature  cover  temperature 
and  the  receiver  temperature  to  the  new  value  calculated. 
The  pressure  drop  through  the  collector  and  the  power 
required  to  pump  the  coolant  through  the  collector  are 
calculated  in  lines  1430  to  1480,  as  specified  in  [Refs.  6 
and  12] .   The  energy  required  to  pump  the  coolant  through 
the  collector  is  then  subtracted  from  the  energy  gained  by 
the  collector  to  calculate  the  available  collector  energy, 
Qa,  line  14  90.   Qa  is  then  summed  for  each  hour  that  the 
analysis  is  performed.   As  the  analysis  is  performed  for 
one  day  each  month,  the  summation  of  the  available  collector 
energy  is  then  multiplied  by  thirty  to  obtain  a  yearly  heat 
gain,  Qy,  line  1590.   The  cost  of  an  equivalent  amount  of 
natural  gas  is  then  calculated.   Using  present  worth  analysis 
as  described  in  [Ref .  13] ,  a  life-cycle  savings  is  then 
calculated  and  the  initial  manufacturing  cost  of  the 
collector  is  calculated  and  subtracted  from  the  life-cycle 
savings.   This  result  is  then  divided  by  the  life-cycle 
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useful  energy  gain  by  the  collector  and  multiplied  by  10 6 
to  obtain  the  life-cycle  fuel  savings  used  as  the  objective 
function  for  the  optimization  program.   The  constraints  on 
the  design  are  then  calculated  in  lines  1660  to  1685.   Lines 
1710  to  1815  contain  printout  specifications  for  the  first 
and  the  last  time  that  the  NISCO  subprogram  is  called  by 
the  optimization  program.   The  values  of  the  objective 
function  and  the  constraints  are  then  passed  along  with 
program  control  back,  to  the  optimization  program. 
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IV.   RESULTS 

A  major  portion  of  the  thesis  work  was  devoted  to  the 
development  of  a  reliable  optimization  program  that  would 
optimize  a  wide  variety  of  problems.   Following  the  accom- 
plishment of  this  goal,  a  subprogram  was  developed  to 
model  a  nonimaging  concentrating  compound  parabolic  trough 
solar  collector. 

A.   RESULTS  OF  THE  DESOP  PROGRAM  DEVELOPMENT 

In  developing  the  DESOP  numerical  optimization  program, 
four  standard  numerical  optimization  test  problems  were 
used.   The  four  test  problems  provided  a  wide  variety  of 
different  numerical  problems  with  which  the  optimization 
program  had  to  deal.   The  four  test  problems  are  listed  in 
Appendix  F.   A  major  goal  of  this  thesis  was  to  optimize 
all  four  problems  using  default  optimizer  control  variables 
The  results  of  the  DESOP  program  optimizing  the  four  test 
problems  are  given  in  Appendix  A.   The  DESOP  program  was 
able  to  make  significant  design  improvements  in  all  four 
test  problems  using  default  optimizer  control  parameters. 
When  the  optimizer  control  parameters  were  adjusted  for 
each  individual  problem,  the  DESOP  program's  performance 
was  improved  in  all  four  cases.   Further  experimenting 
with  the  optimizer  control  parameters  could  lead  to  an 
even  better  performance  of  the  DESOP  program. 
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B.   RESULTS  OF  THE  NISCO  SUBPROGRAM 

A  listing  of  the  NISCO  subprogram  is  included  as 
Appendix  G.   There  are  numerous  comment  statements  in  the 
subprogram.   The  reader  is  encouraged  to  look  closely  at 
the  subprogram.   The  results  of  the  NISCO  subprogram  may 
be  found  in  Appendix  B.   The  first  section  of  Appendix  B 
details  the  design  specifications  that  were  chosen  for  the 
solar  collector  model.   Three  different  daily  heat  capa- 
cities were  specified  10000,  30000  and  50000  Btu,  and  the 
NISCO  subprogram  was  used  to  find  the  optimum  design  for 
each.   The  second  section  of  Appendix  B  gives  the  initial 
design  and  final  designs  for  the  three  solar  collector 
capacities.   In  each  case  the  DESOP  program  was  able  to 
significantly  improve  the  design.   The  instantaneous 
efficiencies  for  the  final  designs  are  within  a  few 
percent  of  the  instantaneous  efficiency  reported  by  Kreith 
and  Kreider  [Ref.  6]  for  a  slightly  different  nonimaging 
concentrating  compound  parabolic  trough  solar  collector 
operating  under  slightly  different  atmospheric  conditions. 
It  is  interesting  to  note  that  the  larger  capacity  solar 
collector  had  the  best  instantaneous  efficiency  and  also 
the  highest  life-cycle  fuel  savings.   In  all  three  cases 
the  optimum  incident  acceptance  angle  was  found  to  be 
18.04  -  18.05  degrees  and  the  optimum  truncation  angle 
was  found  to  lie  within  the  range  of  184.55  to  190.0  degrees 
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Due  to  the  fact  that  the  objective  function  is  weakly 
linked  to  the  design  capacity  of  the  solar  collector,  the 
final  daily  heat  gain  is  somewhat  higher  than  the  minimum 
set.   If  a  stronger  link  were  to  be  established  between 
the  solar  collector  capacity  and  the  objective  function, 
the  collector  design  would  be  driven  closer  to  the  stated 
minimum  daily  average  heat  gain.   Also,  if  the  convergence 
criteria  is  tightened,  the  optimizer  will  take  longer  to 
reach  the  minimum  but  will  reach  a  final  design  where  the 
average  daily  heat  gain  is  closer  to  the  minimum  set. 
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APPENDIX  A 
DESOP  TEST  PROGRAM  RESULTS 
This  appendix  contains  the  results  of  the  four  test 
programs  that  were  used  to  develop  the  DESOP  numerical 
optimization  program.   For  each  design  the  initial  design, 
the  true  optimum  design,  the  DESOP  results  using  default 
control  parameters,  and  the  DESOP  results  using  adjusted 
control  parameters  are  given.   The  four  test  programs  may 
be  found  in  Appendix  F. 
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DESOP  TEST  PROGRAM 


ANALIZ  Suboroaram  :  BANNA 


Initial  Design: 

Desian  Variables: 

X(l)  =  -1.2 

X(2)  =  1.0 
Objective  Function: 

Obj  =10.8 

Opj  =10.8 
Side  Constraints  Violated: 

N/A 
Constraints  Violated: 

N/A 


True  Optimum: 

Desicm  Variables: 

X(l)  =  1.00 

X(2)  =  1.00 
Objective  Function: 

Obj  =  4.00 

Opj  =  4.00 
Side  Constraints  Violated 

N/A 
Constraints  Violated: 

N/A 


DESOP  Results 


Default  Control  Parameters: 
Design  Variables: 

X(l)  =  Q.768 

X(2)  =  0.578 
Objective  Function: 

Obj  =  4.0  5 

Opj  =  4.05 
Side  Constraints  Violated 

N/A 
Constraints  Violated: 

N/A 
#  of  Function  Evaluations 

96 


Adj.usted  Control  Parameters 
Design  Variables: 

X(l)  =  0.791 

X(2)  =  0.606 
Objective  Function: 

Obj  =  4.047 

Opj  =  4.047 
Side  Constraints  Violated 

N/A 
Constraints  Violated: 

N/A 
#  of  Function  Evaluations: 

91 


NOTE:   The  BANNA  subprogram  has  2  design  variables  and  no 
constraints . 
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DESOP    TEST    PROGPAM 


ANALIZ   Subprogram    :    Rcsen-Suzuki 

Initial   Design: 

Design  Variables: 

X(l)    =   1 

X(2)    =    1 

X(3)    =    1 

X(4)  =  1 
Objective  Function: 

Obj  =  31 

Opj  =  31 
Side  Constraints  Violated: 

N/A 
Constraints  Violated: 

None 


True  Optimum: 

Design  Variables: 

X(l)  =  o.o 

X(2)  =  1.0 

X(3)  =  2.0 

X(4)  =  -1.0 
Objective  Function: 
^  .Obj  -  6.00 

Opj  =6.00 
Side  Constraints  Violated 

N/A 
Constraints  Violated 

None 


DESOP  Results: 


Default  Control  Parameters 
Design  Variables : 

X(l)  =  4.72E-02 

X(2)  =  0.998 

X(3)  =  1.98 

X(4)  =  -1.00 
Objective  Function: 

Obj  =  6.088 

Opj  =  6.093 
Side  Constraints  Violated 

N/A 
Constraints  Violated: 

GC3)  =  0.000527 
#  of  Function  Evaluations 
392 


Adjusted  Control  Parameters 
Design  Variables : 

X(l)  =  -5.167E-03 

X(2)  =  1.019 

X(3)  =  1.999 

X(4)  =  -0.9951 
Objective  Function: 

Obj  =  5.9998 

Opj  =  6.007 
Side  Constraints  Violated 

N/A 
Constraints  Violated: . 

G(3)  =  0.00232 
#  of  Function  Evaluations; 
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NOTE:  The  Rosen-Suzuki  subprogram  has  four  design  variables 
and  three  constraints . 
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DESOP  TEST  PROGRAM 


ANALIZ  Subcroaram  :  T7VAR 

Initial  Design: 

Design  Variables: 


True  Optimum: 

Design  Variables 


X(l)  =  1 

X(l)  =  3 

X(2)  =  1 

X(2)  =  0 

X(3)  =  1 

X(3)  =  0 

X(4)  =  1 

X(4)  =  1 

X(5)  =  1 

X(5)  =  0 

X(6)  =  1 

X(6)  =  0 

X(  )  =  1 

X(7)  =  0 

Side  Constraints  Violated: 

Side  Constraints  Violated 

None 

None 

Objective  Function: 

Objective  Function 

Obj  =  -203 

Obj  =  -190 

Opj  =  1007 

Opj  =  -190 

Constraints  Violated: 

Constraints  Violated: 

G(l)  =  11 

None 

DESOP  Results  : 

Default  Control  Parameters 
Design  Variables 


Adjusted  Control  Parameters 
Design  Variables: 


X(l)  =  1.42 

X(l)  =  2.69 

X(2)  =  0.515 

X(2)  =  -7.60E-04 

X(3)  =  0.376 

X(3)  =  1.04E-02 

X(4)  =  0. 844 

X(4)  =  1.82 

X(5)  =  2.13E-02 

X(5)  =  2.89E-03 

X(6)  =  9.84E-03 

X(6)  =  2.91E-03 

X(7)  =  0.638 

X(7)  =  0.113 

Side  Constraints  Viola 

.ted: 

Side  Constraints  Violated 

None 

1 

Objective  Function: 

Objective  Function: 

Obj  =  -142 

Obj  =  -178.6 

Opj  =  -142 

Opj  =  -179.1 

Constraints  Violated: 

Constraints  Violated: 

None 

None 

#  of  Function  Evaluations : 

#.of  Eunction  Evaluations 
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808 

NOTE:.  The  T7VAR  subprogram  ha 


3  seven  design  variables 


ana  one 


constraint . 


APPENDIX   B 


NISCO  DESIGN  RESULTS 


DESIGN  SPECIFICATIONS 

A.  GEOGRAPHIC 

Location 

Solar  position  and 
intensity  levels. 

Wall  azimuth  angle 
Cloud  Cover 
Vapor  Pressure 
Tamb 

B.  COLLECTOR  MATERIALS 

RECEIVER 

Absoptivity 
Thermal  Emissivity 
Reflectivity 

REFLECTOR 

Reflectivity 


40  degrees  North  Lattitude 
ASHRAE  Handbook  of  Fundamentals 

standard  solar  radiation  design, 
Ref. 
0  deg. 
0  % 

3  mm  Hg 
ASHRAE  Handbook  of  Fundamentals 

daily  norms  for  San  Francisco,  Ca 

Ref. 


Oxidized  Copper 

0.93 

0.40 

0.07 

Vacuum  deposited  Aluminum  on  resin 
0.89 


COVER  Double  strength  window  glass 

Specular  Absorptance   ASHRAE  standards  for  D.S.  window 
Specular  Transmisivity   glass.   Ref. 


Absorptance  av. 
Transmittance  av. 
Reflectivity  av. 
Thermal  Emissivity 

COOLANT 

Specific  gravity 
Specific  heatt 
Viscosityt 
Inlet  Temperature 


0.03 
0.80 
0.17 
0.94 

Therminol  55 

0.87 

4.9E-4*Tr+0.4036   Btu/(lbm  F) 
6.71955E-4*(-0.053*Tr+32.3)  lbm/(ft  sec) 
100  deg.  F 


t  Values  given  as  a  function  of  average  receiver  temperature, 
Tr,  in  degrees  F. 
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C.   ECONOMIC 

Life-Cycle  20  years 

Fuel  Cost  (Natural  Gas)  8.14   $/ 10   Btu 

Annual  Fuel  Inflation  Rate  0.11 

Monetary  Inflation  Rate  0.10 

Collector  Cost  22.5   $/ft3   Ref.  4 
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NISCO  DESIGN   -   10000  Btu/DAY  SOLAR  COLLECTOR 


INITIAL  DESIGN 

Design  Variables: 

Incident  acceptance  half  angle 

Truncation  angle 

Receiver  radius 

Collector  length 

Coolant  mass  flow  rate 
Design  Features : 

Concentration  ratio 

Collector  aperature  area 

Collector  depth 

Coolant  velocity 

Average  daily  heat  gain 

Maximum  coolant  temperature 

Instantaneous  collector  efficiency 

Initial  cost 

Life-cycle  cost  savings 


25.00  deg 

230.0  deg 

0.41  in 

20.00  ft 

10  8.7  lbm/hr 

2.20 

9.56  ft2 

4.12  in 

0.15  ft/sec 

11,400  Btu 

137  deg  F 
0.70 
$203.00 

$7.78  /106  Btu 


FINAL  DESIGN 

Design  Variables : 

Incident  acceptance  half  angle 

Truncation  angle 

Receiver  radius 

Collector  length 

Coolant  mass  flow  rate 
Design  Features : 

Concentration  ratio 

Collector  aperature  area 

Collector  depth 

Coolant  velocity 

Average  daily  heat  gain 

Maximum  coolant  temperature 

Instantaneous  collector  efficiency 

Initial  cost 

Life-cycle  cost  savings 


18.05  deg 
190.0   deg 
0.52  in 

18.92  ft 
124.3   lbm/hr 

1.65 
8.48  ft2 
1.02  in 
0.11  ft/sec 
10,600    Btu 
131    deg  F 
0.74 
$159.00 

$8.15  /106  Btu 
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NISCO  DESIGN  -  30000  Btu/DAY  SOLAR  COLLECTOR 


INITIAL  DESIGN: 

Design  variables: 

Incident  acceptance  half  angle 

Truncation  angle 

Receiver  radius 

Collector  length 

Coolant  mass  flow  rate 
Design  features: 

Concentration  ratio 

Collector  aperature  area 

Collector  depth 

Max.  coolant  temperature 

Instantaneous  collector  efficiency 

Average  daily  heat  gain 

Initial  cost 

Life  cycle  savings 


15.00 

deg. 

230.0 

deg. 

1.00 

in. 

50.0 

ft. 

100.0 

lbm/hr 

2.95 

77.35 

ft2 

12.62 

in. 

366.59 

F 

0.64 

83,400 

Btu 

$1,628.90 

/106    B 

$7.54 

Btu 


FINAL  DESIGN: 

Design  variables: 

Incident  acceptance  half  angle 

Truncation  angle 

Receiver  radius 

Collector  length 

Coolant  mass  flow  rate 
Design  features: 

Concentration  ratio 

Collector  aperature  area 

Collector  depth 

Coolant  velocity 

Max.  coolant  temperature 

Instantaneous  collector  efficiency 

Average  daily  heat  gain 

Initial  cost 

Life-cycle  savings 


18.04 

deg. 

184.55 

deg. 

1.28 

in. 

26.43 

ft. 

425.23 

lbm/hr 

1.54 

27.35 

ft? 

1.78 

in. 

6.04E-02 

ft/sec 

129.4 

F 

0.75 

34,473 

Btu 

$509.98            , 
$8.19    /10D 

Btu 


Run  time   -   approx.    5  hours 
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APPENDIX  C 
DESOP  USER'S  MANUAL 

A.   INTRODUCTION 

The  numerical  optimization  package  of  programs,  currently 
consisting  of  OPCON  and  DESOP,  provide  the  capability  for 
finding  the  optimum  design  of  a  system  mathematically  modeled 
using  multiple  variables,  on  the  Hewlett-Packard  9845A  desk- 
top computer.   OPCON  is  an  executive  program  which  provides 
the  user  with  the  following:   a  primary  point  of  contact  with 
the  computer  from  which  to  access  the  optimization  program 
DESOP,  a  standard  formatted  input  for  design  variables,  side 
constraints  on  the  design  variables,  optimizer  control  varia- 
bles, and  organization  of  the  optimization  process.   DESOP 
is  a  numerical  optimization  program,  which  when  coupled  to 
a  user  supplied  design  analysis  program,  will  optimize  the 
design.   DESOP  will  allow  the  user  to  monitor  the  optimiza- 
tion process  as  it  is  taking  place.   After  monitoring  the 
optimization  process,  the  user  may  choose  to  change  the 
optimizer  control  variables  and/or  his  design  starting 
point  to  more  efficiently  or  more  accurately  reach  the 
design  optimum.   After  an  optimization  has  been  performed, 
DESOP  will  reload  the  OPCON  program  and  return  the  user  to 
the  OPCON  program.   The  OPCON  program  will  then  offer  the 
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user  the  following  three  choices:   to  optimize  the  same 
design,  optimize  a  different  design  or  terminate  the 
program. 

Figure  11  illustrates  the  overlay  method  used  in 
loading  the  OPCON  and  DESOP  programs. 
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PROGRAM  OVERLAY  STRUCTURE 


MASS  STORAGE 


OPCON 


Input  information 
to  files 


DESOP 


Read  information 
from  files 


Optimization  Program 


User' s 

Design  Analysis 

Subprogram 


/ 


/ 


/ 


/ 


It. 


OPCON 


DESOP 


ANALYSIS  #1 
Subprogram 


ANALYSIS  #1 
Data  Files 


ANALYSIS  #2 
Subprogram 


ANALYSIS  #2 
Data  Files 


ANALYSIS  #N 
Subprogram 


ANALYSIS  #N 
Data  Files 


Figure  11.   Basic  Program  Relationships 
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B.   THE  CONCEPT  OF  NUMERICAL  OPTIMIZATION 

Consider  the  problem  of  designing  the  cantilevered  beam 
shown  in  Figure  1.  The  design  task  may  be  broken  down  into 
three  major  parts.  First,  the  objective  of  the  design  must 
be  determined,  which  in  this  case  is  to  minimize  the  weight 
of  the  beam  required  to  support  the  concentrated  tip  load  P. 
Second,  any  physical  constraints  that  may  effect  the  design 
must  be  determined.  Thirdly,  any  limits  which  exist  on  the 
design  variables  must  be  determined.  The  design  problem 
may  then  be  reduced  to  a  system  of  equations  as  follows: 

Minimize  the  volume  (V) 
V  =  B-H-L 

Subject  to: 

Bending  stress  (a,  ) 

o 

6  •  P«  L 
°b  =  V^  -  2000°  psi 

Shear  stress  ('  ) 


v  =  Ath  -   1000°  ?si 


Deflection  under  load  (:) 

4-P-L3    ,  • 
5  =  —  _  ..  ;  _  1  inch 

r.  •  ,  -.  •  h  ' 


Ci  •  O 
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P=10000 


h— I 


Figure  1.   Cantilevered  Beam  Design  Problem 
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With  geometric  constraints  such  that: 
0.5  ^  B  <  0.5 
1.0  <    H  <  20.0 
H/B  <  10.0 
At  this  time  the  following  definitions  are  introduced: 

Objective  Function  -  The  parameter  which  is  to  be 
minimized  or  maximized  during  optimization.   The  objective 
function  always  occurs  on  the  left  side  of  the  equation 
unless  it  is  also  used  as  a  design  variable.   An  objective 
function  may  be  either  linear  or  nonlinear,  implicit  or 
explicit,  but  must  be  a  function  of  the  design  variables. 

Design  Variable  -  Any  parameter  which  the  optimiza- 
tion process  is  allowed  to  change  in  order  to  improve  the 
design.   Design  variables  appear  only  on  the  right  hand 
side  of  equations  in  the  analysis  program. 

Inequality  Constraint  -  Any  parameter  which  must 
not  exceed  specified  bounds  for  the  design  to  be  acceptable. 
Constraint  functions  always  appear  on  the  left  side  of 
equations.   A  constraint  may  be  linear,  nonlinear,  implicit 
or  explicit,  but  must  be  a  function  of  one  or  more  design 
variables. 

Equality  Constraint  -  Any  parameter  which  must  equal 
a  specified  value  for  the  design  to  be  acceptable.   The 
same  rules  apply  to  equality  constraints  as  inequality 
constraints . 
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Side  Constraint  -  Any  upper  or  lower  bound  placed 
upon  a  design  variable.   Side  constraints  are  usually  not 
included  in  the  system  of  equations  that  comprise  an 
analysis  program.   Instead  they  are  usually  included  as 
part  or  the  data  input  to  the  optimization  program. 

Analysis  Code  -  The  system  of  equations  utilizing 
the  design  variables  which  are  used  to  calculate  the  objective 
function  and  the  constraints  of  a  particular  design  problem. 

The  general  optimization  problem  may  thus  be  stated 
mathematically  as: 

Find  the  set  of  design  variables  X.  where  i  =  l,2,...,n 
which  will: 

Minimize  the  objective  function  (Obj) 

Obj  =  f(X) 
Subject  to: 

Inequality  constraints  (G) 

G.  (X)  <  0       j  =  1,2, ...  ,m 

Equality  constraints  (H) 

H. (X)  =  0        j  -  1,2,... ,1 

Side  constraints 

XT  <    X.  £  Xu     i-=  1,2,. . .  ,n 
11     i 

Returning  to  the  cantilevered  beam  problem,  it  may  be 

stated  in  the  standard  format  as  follows: 

Le-t  X(l)  =  3,  X(2)  =  H,  and  Obj  =  Vol  =  B-H-L 

Then  minimize  Obj  =  Vol 
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Subject  to: 

G(1)  =  20M0  "  X  -  ° 

G(2)  =  loWO  '  l    ~    ° 
G(3)=o-l<0 

G(4)  =  |  -  10  <    0 

With  side  constraints: 

X(l)  -1  =  0.5 

X(1)U  =  5.0 

X(2)1  =  1.0 

X(2)U  =  20.0 
It  is  thus  fairly  simple  and  straightforward  to  perform 
an  analysis  on  a  particular  beam  for  a  given  B  and  H. 
Successive  analyses  may  be  performed  on  the  cantilevered 
beam  by  solving  the  above  system  of  equations.   It  is  de- 
sirable to  automate  the  successive  solutions  and  to  direct 
the  solutions  such  that  each  solution  is  a  better  design 
than  the  last.   One  approach  for  doing  so,  and  the  one 
used  by  DESOP  is  to  proceed  as  follows:   Start  with  initial 
values  for  B  and  H.   Solve  the  above  set  of  equations  to 
find  the  objective  function  Obj  and  to  see  if  any  constraints 
are  violated.   A  pseudo  objective  function  is  created  to 
represent  designs  when  constraints  are  violated.   If  a 
constraint  is  violated*,  3.   penalty  is  added  to  Obj  to  form 
a  penalized  objective  function  Op j .   The  gradient  of  the 
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penalized  objective  function  at  the  initial  design  may  be 
found  by  taking  the  first  partial  derivative  of  Opj  with 
respect  to  the  design  variables.   The  gradient  of  the 
penalized  objective  function  defines  the  direction  of 
steepest  ascent.   In  the  case  of  the  cantilevered  beam,  it 
is  desired  to  minimize  the  objective  function;  therefore, 
the  greatest  improvement  in  design  may  be  achieved  by  moving 
in  the  negative  gradient,  or  steepest  descent  direction. 
From  the  initial  design  point  a  search  is  performed  in  the 
steepest  descent  direction  for  the  minimum  value  of  Opj  in 
that  direction.   At  the  new  minimum,  the  gradient  of  the 
penalized  objective  function  is  again  determined  and  a 
search  is  performed  in  a  conjugate  direction  until  a  second 
minimum  is  found.   Successive  iterations  are  performed 
until  the  gradient  is  found  to  be  zero  or  each  successive 
iteration  produces  a  sufficiently  small  change  in  Opj  such 
that  for  all  practical  purposes  the  minimum  has  been  found. 
At  this  time  the  penalty  function  is  increased.   If  the 
design  is  in  a  region  where  there  are  no  constraints 
violated  an  increase  in  the  penalty  function  will  not 
change  the  value  of  Opj.   If  on  the  other  hand  the  design 
is  in  an  infeasible  region  where  there  are  one  or  more 
constraints  violated,  Opj  will  be  increased,  and  the 
search  for  a  new  minimum  will  commence.   If  the  minimum 
of  the  objective  function  exists  in  the  infeasible  region, 
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the  minimum  value  for  the  objective  function  in  the  feasible 
region  will  be  approached  from  the  infeasible  region  as  the 
penalty  function  is  increased.   The  design  improvement 
process  will  terminate " when  a  zero  gradient  is  found  or 
successive  iterations  produce  a  sufficiently  small  change 
in  the  value  of  Opj  and  an  increase  in  the  penalty  function 
causes  no  change  in  Opj . 

The  minimum  thus  reached  by  the  optimization  process  is 
a  minimum  with  respect  to  the  penalized  objective  surface 
immediately  surrounding  the  final  design  point.   The  opti- 
mization process  cannot  distinguish  between  local  and  global 
minimum  points.   It  is  thus  good  engineering  practice  to 
run  several  optimizations  for  a  particular  design  problem 
from  several  different  initial  design  points.   If  optimiza- 
tions performed  from  different  initial  design  points  converge 
on  the  same  minimum  point,  that  point  is  probably  a  global 
minimum.   If  on  the  other  hand  two  or  more  minimums  are 
found,  there  may  be  local  minimums  located  in  the  design 
space  being  considered  and  care  must  be  taken  to  find  the 
global  minimum. 

C   THE  DESOP  NUMERICAL  OPTIMIZATION  PROGRAM 

The  DEsktop  Sequential  unconstrained  minimization 
technique  Optimization  Program  was  developed  using  the 
basic  optimization  approach  outlined  in  Section  II. A. 
A  copy  of  the  program  is  included  in  Appendix  E.   The  major 
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program  structure  is  shown  in  Figure  2.   The  following 
discussion  will  refer  to  Figure  2  and  describe  the  major 
features  of  the  program. 

1.   Basic  Program  Execution 

The  DESOP  program  begins  execution  when  it  is 
loaded,  linked  to  the  user's  analysis  subprogram,  and  the 
program  is  instructed  to  run  by  the  OPCON  program.   The 
above  actions  are  automatically  performed  by  the  OPCON 
program.   DESOP  is  loaded  into  the  computer  by  an  overlay 
process.   Therefore  no  variables  can  be  directly  transferred 
between  the  DESOP  and  OPCON  programs.   DESOP  begins  execu- 
tion by  reading  the  optimizer  control  variables  and  the 
design  variables  that  were  input  using  the  OPCON  program 
and  saved  to  a  mass  storage  device.   The  program  then  sets 
Icalc  equal  to  one  and  evaluates  the  objective  function 
and  constraints  at  the  initial  design  point.   Icalc  is  a 
flag  provided  the  user  to  key  user  specified  output  on  the 
initial  and  final  design  analysis.   DESOP  will  provide  the 
user  with  a  hard  copy  output  of  the  initial  design  variables, 
the  value  of  the  constraints,  the  objective  function  and 
the  penalized  objective  function.   The  user  then  has  the 
option  of  continuing  with  the  DESOP  program  to  optimize 
his  analysis  subprogram  or  to  return  to  the  OPCON  program 
to  change  one  or  more  of  the  input  parameters. 
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INPUT  X0 


EVALUATE  Obj  &  G 


FIND  THE  GRADIENTS 
OF  Obj  AT  X 


CALCULATE  THE 
STEEPEST  DESCENT 
SEARCH  DIRECTION 


CALCULATE  THE 
FLETCHER  REEVES 
SEARCH  DIRECTION 


ESTIMATE  a  FOR  MIN.  Obj 
IN  1-D  SEARCH 


FIND  a  FOR  MIN.  Obj 
BY  EXPANDING  SEARCH, 
METHOD  OF  GOLDEN  SECTION, 
&  POLYNOMIAL  INTERPOLATION 


MOVE  TO  NEW  X 


CHECK  CONVERGENCE  OF  A Obj 


INCREASE  PENALTY  FUNCTION 
RE CHECK  CONVERGENCE 


IF  MET 


END 


Figure  2.   DESOP  Flow. Diagram 
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Proceeding  with  the  optimization,  there  are  two 
major  loops  in  the  optimization  program.   The  outer  loop 
increases  the  penalty  function  when  the  inner  loop's  con- 
vergence criteria  have  been  met.   A  convergence  test  is 
then  performed  by  the  outer  loop.   If  the  convergence 
criteria  is  met,  the  optimization  process  is  considered 
finished.   If  the  convergence  criteria  for  the  outer  loop 
is  not  met,  program  execution  is  returned  to  the  inner 
loop.   The  inner  loop  performs  successive  iterations 
searching  for  the  minimum  of  the  penalized  objective 
function  with  no  increase  in  the  penalty  function  taking 
place.   When  the  inner  loop's  convergence  criteria  have 
been  met  program  execution  is  transferred  to  the  outer  loop. 

Execution  of  the  program  while  in  the  inner  loop 
proceeds  as  follows:   First,  the  gradient  of  the  penalized 
objective  function  is  calculated  by  subroutine  GRAD.   The 
program  then  computes  a  search  direction  using  either  a 
steepest  descent  method  or  the  method  of  conjugate  directions 
developed  by  Fletcher  and  Reeves  [Ref .  1] .   Once  a  search 
direction  is  established  the  optimizer  attempts  to  locate 
the  minimum  value  of  the  penalized  objective  function  in 
the  search  direction.   This  process  if  referred  to  as  the 
one-dimensional  search  and  is  illustrated  in  Figure  3. 
The  efficiency  and  accuracy  to  which  the  one-dimensional 
search  for  the  minimum  of  the  penalized  objective  function 
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is  accomplished,  is  the  key  to  successful  sequential  un- 
constrained minimization  technique  numerical  optimization. 
The  one-dimensional  problem  may  be  expressed  in  terms  of 
the  penalized  objective  function,  Opj ,  and  the  amount  of 
movement,  .,  in  the  search  direction,  S.   First  the  slope 
of  Opj  with  respect  to  «  is  calculated.   An  initial  "guess" 
of  how  far  to  move  is  made  using  subroutine  ALPGES.   The  a 
which  corresponds  to  the  minimum  value  of  Opj  in  the  one- 
dimensional  search  is  then  calculated  using  subroutine 
ALPBND  and  subroutine  QUEBIC.   The  minimum  value  of  Opj 
thus  found  is  then  compared  to  the  previous  value  of  Opj 
for  convergence  using  subroutine  CONVRG.   If  convergence 
is  not  met,  execution  returns  to  the  start  of  the  inner 
loop.   If  convergence  is  met,  execution  returns  to  the 
outer  loop. 

When  the  convergence  criteria  have  been  met  for 
both  the  inner  and  outer  loops,  the  program  proceeds  to 
set  ICALC  to  three  as  a  flag  for  user  generated  output  for 
the  final  design.   DESOP  then  provides  the  user  with  a 
hard  copy  output  of  the  final  design  variables,  objective 
function  value,  penalized  objective  value,  constraint 
values,  the  number  of  inner  loop  iterations,  the  number  of 
times  the  analysis  subprogram  was  called  and  the  final  value 
of  the  penalty  function.   The  OPCCN  program  is  then  overlayed 
over  the  CESOP  program  and  program  execution  is  returned  to 
the  OPCON  program. 
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2.   Finding  the  Search  Direction 

The  first  step  in  finding  the  search  direction,  S, 

is  to  determine  the  slope  of  Opj  at  the  present  design 
point.   The  forward  finite  difference  method  is  used  where: 

3P       F(xi  +  AV  -  F(V   q 

dX.  LX.  i 

1  1 

i  =  1,2,... Ndv 

As  3F/jX.  gives  the  direction  of  positive  slope,  the  search 
direction  is  the  negative  of  3F/3X. .   The  first  search  is 
performed  using  the  steepest  descent  as  found  above  using 
the  following  relation: 

X  r  =  X .  +  aS  . 
11      l 

where  alpha  is  the  distance  moved  in  the  S  direction.   When 
a  minimum  is  obtained  along  the  direction  of  steepest 
descent,  a  new  Fletcher-Reeves  conjugate  search  direction 
[Ref.  1]  is  calculated  at  the  new  Design  point  using  the 
following  relations: 

s:  =  J£-  +  bs. 

1    a  X .       1 


Ndv 


3F'  i  2 


B  =  1=1      X 


Nav   r_2£_i 
Z    L  X  .  J 
i=l 


where  the  prime  denotes  values  for  the  present  iteration 
and  the  non-prime  variables  indicate  values  for  the  previous 
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iteration.   A  one-dimensional  search  is  then  performed  in 
the  new  search  direction.   Searches  are  continued  using 
the  conjugate  direction  method  for  Ndv  +  1  iterations, 
where  Ndv  is  the  number  of  design  variables.   The  search 
process  is  then  restarted  using  the  steepest  descent  method. 
The  reason  for  incorporating  the  conjugate  search  method 
is  that  the  steepest  descent  method  when  traversing  a  design 
surface  with  a  curved  valley  will  tend  to  zigzag  from  one 
side  of  the  design  surface  valley  to  the  other  making  very 
little  progress  as  is  illustrated  in  Figure  4.   The  conju- 
gate direction  method  is  much  more  efficient  in  traversing 
such  a  design  surface.   However,  as  the  conjugate  direction 
method  is  additive  upon  previous  searches,  it  has  a  tendency 
to  decrease  in  effectiveness  with  each  successive  search 
owing  to  the  accumulation  of  numerical  "noise."   For  that 
reason  the  search  process  is  restarted  with  the  steepest 
descent  method  every  Ndv  +  1  iterations,  or  when  the  con- 
jugate direction  predicts  a  positive  slope.   The  search 
direction  is  normalized  to  avoid  inaccuracies  caused  by 
numerical  ill-conditioning. 

3.   Estimating  an  Initial  Value  for  Alpha 

The  initial  estimate  for  alpha  is  made  in  the 
following  manner:   First,  the  slope  of  Opj  in  the  search 
direction  is  calculated  as  the  sum  of  each  of  the  products 
of  the  gradients  times  the  search  direction.   Then  the 
slope  of  Opj  in  the  search  direction  is  divided  by  the  value 
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of  Opj .   This  value  is  then  multiplied  by  an  improvement 
percentage  in  Opj .   This  first  estimate  is  then  applied  to 
a  series  of  conditional  tests  to  determine  the -validity  of 
the  estimate  with  respect  to  the  slope  of  Opj  and  the  magni- 
tude of  the  design  variables.   Lastly,  the  estimate  for 
alpha  is  checkec  to  see  if  it  violates  any  side  constraints. 
If  it  does,  the  value  of  the  estimate  for  alpha  is  reduced 
until  the  side  constraints  are  no  longer  violated. 

4.   Calculating  Alpha 

The  calculation  of  alpha  is  the  most  critical  al- 
gorithm in  the  DESOP  program  in  providing  reliable  optimizer 
operation.   The  ability  to  accurately  and  efficiently  find 
the  minimum  of  the  penalized  objective  function  in  the  one- 
dimensional  search  affects  directly  the  operation  of  the 
optimizer.   Figure  5  illustrates  the  zigzag  phenomenon 
which  occurs  if  alpha  is  not  accurately  found.   The  zigzag 
phenomenon  is  caused  by  the  fact  that  the  optimizer  in  per- 
forming the  forward  finite  difference  for  calculating  the 
search  direction  perturbs  the  design  vector  a  very  small 
amount.   As  such  the  optimizer  can  only  "see"  the  design 
surface  that  is  immediately  adjacent  to  the  design  point. 
Therefore,  if  the  minimum  is  not  found  in  the  one-dimensional 
search,  the  optimizer  will  converge  very  slowly  on  the 
minimum. 

There  are  tew  major  sections  to  the  ALPBND  subroutine. 
The  first  section  attempts  to  find  the  minimum  value  of  Opj 

86 


Figure  5.   The  Zigzag  Phenomenon 
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using  an  expanding  search  technique.   The  first  move  is 
the  amount  predicted  by  the  ALPGES  subroutine.   If  the 
minimum  is  not  bracketed  by  the  first  move,  additional 
moves  are  made.   Each  move  is  larger  than  the  last.   The 
size  o£  the  move  is  increased  each  time  by  an  amount  equal 
to  the  size  of  the  last  move  divided  by  the  lower  Golden 
Section  fraction,  where  the  Golden  Section  fractions  are: 

Fl 2 

•  F   =  *'!  -  l 
r2      2 

The  lower  Golden  Section  fraction,  F.  ,  is  used  so  that  the 
interval  will  be  consistent  with  the  Golden  Section  search 
in  the  second  section  of  the  ALPBND  subroutine.   The  expand- 
ing search  is  continued  until  the  minimum  value  of  the  ob- 
jective function  has  been  bracketed. 

Once  the  minimum  is  bracketed,  a  Golden  Section 
search  is  performed  to  reduce  the  bracketing  interval  on 
the  minimum  by  an  amount  such  that  when  the  two  end  points 
of  the  interval  are  taken  with  two  points  internal  to  the 
interval  and  a  cubic  is  passed  through  the  four  points,  the 
cubic  will  accurately  predict  the  minimum  of  the  penalized 
objective  function.   Himmelblau  in  [Ref.  2]  states  that  the 
Golden  Section  search  method  of  reducing  the  interval 
around  the  minimum  of  Opj  is  the  most  effective  of  the 
reducing  techniques  studied.   Golden  Section  search  is  based 
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on  the  splitting  of  a  line  into  two  segments  known  in 
ancient  times  as  the  "Golden  Section. "   The  ratio  of  the 
whole  line  to  the  larger  segment  is  the  same  as  the  ratio 
of  the  larger  segment  to  the  smaller  segment.   The  two 
Golden* Section  fractions  are  employed  to  split  the  interval 
bracketing  the  minimum  as  shown  in  Figure  3.   Once  the 
interval  has  been  split,  the  two  values  of  Opj  corresponding 
to  the  internal  points  are  compared  to  find  the  larger  of 
the  two.   The  internal  point  with  the  larger  value  of  Opj 
will  become  the  new  end  point  for  the  interval,  the  remain- 
ing interior  point  will  by  the  fact  that  it  was  determined 
by  a  Golden  Section  fraction,  be  equal  to  the  point  deter- 
mined by  the  other  Golden  Section  fraction.   Thus,  only  one 
new  point  must  be  calculated  to  continue  the  Golden  Section 
search.   The  search  is  continued  in  this  manner  until  the 
vertical  separation  of  the  two  end  points  with  respect  to 
the  interior  points  is  less  than  one  percent.   The  four 
values  of  the  penalized  objective  function  corresponding  to 
the  four  Golden  Section  search  points  are  then  sent  to  a 
cubic  interpolator.   The  cubic  interpolator  will  return  a 
value  for  alpha  to  predict  the  minimum  of  the  penalized 
objective  function,  and  the  minimum  of  the  cubic  function 
that  the  interpolator  has  created.   The  subroutine  ALPBND 
will  then  test  the  predicted  minimum  with  the  minimum  found 
at  the  predicted  alpha.   If  there  is  less  than  a  tenth  of 
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one  percent  difference  between  the  two  values  of  the  objec- 
tive function,  the  point  predicted  by  the  cubic  interpolator 
will  be  accepted  as  the  minimum  and  program  execution  will 
return  to  the  main  program.   If  the  predicted  minimum  is  not 
sufficiently  close  to  the  minimum  at  the  predicted  alpha, 
another  Golden  Section  search  will  be  performed  to  reduce 
the  interval  and  better  localize  the  minimum.   The  four 
points  from  the  reduced  interval  will  then  be  sent  to  the 
cubic  interpolation  subroutine.   This  process  will  continue 
until  either  the  test  for  the  minimum  is  positive  or  the 
interval  has  been  reduced  to  less  than  1E-12.   Program  exe- 
cution will  then  return  to  the  main  program. 
5.   Subroutine  QUEBIC 

Subroutine  QUEBIC  is  used  to  estimate  the  alpha  at 
which  Opj  is  a  minimum  based  on  four  point  cubic  interpola- 
tion.  If  the  function  more  closely  resembles  a  quadratic 
than  a  cubic,  a  three  point  quadratic  interpolation  is 
performed  using  the  three  points  which  bracket  the  minimum. 
If  the  predicted  minimum  is  outside  the  interval  spanned 
by  the  two  end  points  again  a  quadratic  interpolation  is 
performed.   If  the  minimum  still  lies  outside  the  two  end 
points,  the  analysis  returns  to  subroutine  ALPBND,  the 
inverval  bracketing  the  minimum  is  reduced,  and  program 
execution  returns  to  QUEBIC. 
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6.  Convergence  of  the  Penalized  Objective  Function 
The  penalized  objective  function  is  tested  for  con- 
vergence at  the  end  of  each  inner  loop  and  again  at  the  end 
of  the  outer  loop  in  the  main  program.   Convergence  is 
tested  by  calling  subroutine  CONVRG.   There  are  two  criteria 
used  for  testing  for  convergence.   The  first  tests  the  rela- 
tive difference  of  the  value  of  Opj  from  the  present  iteration 
with  the  value  of  Opj  from  the  last  iteration.   The  second 
method  tests  the  absolute  difference  of  the  two  values.   The 
second  method  is  employed  for  cases  when  the  value  of  the 
penalized  objective  function  approaches  zero.   When  conver- 
gence has  been  met  on  two  successive  iterations,  the  penalty 
function  is  increased  by  an  amount  specified  by  the  user  in 
the  executive  OPCON  program.   The  penalized  objective  function 
is  again  tested  for  convergence.   If  convergence  is  still 
met,  the  optimizer  considers  the  present  value  of  the  pena- 
lized objective  function  to  be  a  minimum,  noting  again  that 
numerical  optimization  programs  cannot  differentiate  between 
local  and  global  minimums. 

7.  The  Penalty  Function 

The  purpose  of  the  penalty  function  is  to  increase 
the  value  of  the  objective  function  when  the  design  is  in 
an  infeasible  region.   The  infeasible  region  is  that  region 
where  one  or  more  design  constraints  are  violated.   When  a 
constraint  is  violated,  the  value  of  the  particular  constraint, 


G-,  is  positive.   The  objective  function  is  then  penalized 
as  follows: 


Opj  =  Obj  +  R.G.E 
whe  re : 

R  -  a  multiplication  constant 

E  -  an  exponent  constant 
This  type  of  penalty  function,  one  where  the  penalty  is 
applied  after  the  design  leaves  the  feasible  region,  is 
known  as  an  exterior  penalty  function.   The  exterior  type 
of  penalty  function  was  chosen  over  other  types,  such  as 
the  interior  or  extended  interior  penalty  function.   If  a 
function  is  discontinuous  within  the  design  space  being 
studied,  numerical  difficulties  may  be  encountered  which 
make  performing  an  optimization  of  the  design  difficult. 
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D.   USE  OF  THE  DESOP  NUMERICAL  OPTIMIZATION  PROGRAM 

Once  the  analysis  subprogram  has  been  written  and  SAVED, 
the  user  may  perform  a  numerical  optimization  of  the  design 
analysis.   Figure  11  shows  the  basic  relationships  between 
the  OPCON  and  DESOP  programs  and  the  user  supplied  ANALIZ 
subprograms.   Note  that  the  DESOP  program  overlays  the  OPCON 
program  after  all  data  has  been  input.   To  begin,  LOAD  OPCON 
and  press  run.   Through  a  series  of  self-explanatory  menus, 
the  user  will  be  prompted  to  input  control  variables,  design 
variables  and  execute  the  DESOP  program.   The  following 
menus  appear  in  the  OPCON  program. 

1.  INTRODUCTION  -  A  brief  introduction  stating  the 
purpose  of  the  OPCON  program. 

2.  NEW  OR  EXISTING  ANALYSIS  PROGRAM  -  If  this  is  the 
first  time  that  a  particular  analysis  subprogram  has  been 
run  on  the  tape  or  disk  being  used,  or  if  the  files  from 
a  subsequent  run  have  been  deleted,  the  user  must  enter 
the  response  for  a  new  analysis  subprogram.   For  existing 
programs,  OPCON  will  read  data  from  the  existing  data 
file.   The  data  read  will  be  that  from  the  subsequent  run 
of  the  particular  analysis  subprogram. 

3.  NAME  THE  ANALIZ  SUBPROGRAM  -  The  user  will  input 
the  name  under  which  the  analysis  subprogram  he  wishes  to 
use  was  saved.   If  this  is  a  new  analysis  subprogram 
OPCON  will  create  files  for  saving  the  data  input  during 
the  execution  of  OPCON. 
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4.  OPTIMIZER  SELECTION  -  At  the  present  time  there  is 
only  one  optimization  program  available.   It  is  hoped  that 
at  some  future  date  additional  numerical  optimization 
programs  will  be  added  to  the  optimization  package.   OPCON 
has  been  developed  to  interface  with  multiple  numerical 
optimization  programs.   Details  for  linking  other  numerical 
optimization  programs  to  OPCON  are  included  as  comments  in 
the  OPCON  listing  which  may  be  found  in  Appendix  D. 

5.  INPUT  NDV  AND  NCON  -  For  a  new  analysis  subprogram 
the  user  will  be  asked  to  input  the  number  of  design  varia- 
bles, NDV,  and  the  number  of  constraints  in  the  analysis 
subprogram,  NCON. 

6.  INPUT  CHECK  OF  COMMON  CONTROL  VARIABLES  -  A  menu  of 
the  optimizer  control  parameters  common  to  DESOP,  Feasible 
Direction  type  and  other  future  optimization  programs  is 
displayed.   The  menu  displays  the  variable  name,  its 
minimum,  maximum,  default  and  present  value.   The  variables 
displayed  are: 

NDV  -  The  number  of  independent  design  variables 
used  in  the  analysis  code. 

NCON  -  The  number  of  constraints  in  the  design 
analysis  subprogram. 

IPRINT  -  A  print  control  used  in  the  optimization 
program  to  display  intermediate  results.  Positive  values 
entered  will  print  on  the  CRT.   Negative  values  entered 
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will  print  on  the  thermal  printer.   If  a  zero  is  entered 
there  will  not  be  a  hard  copy  output  of  the  initial  and 
final  results  of  the  optimization  program.   IPRINT  =  0  is 
to  be  used  when  debugging  an  analysis  subprogram  to 
conserve  thermal  paper. 

+1  -  Print  initial  and  final  optimization 
information. 

+2  -  Print  above  plus  the  objective 
function  and  penalized  objective  function  on  each  iteration. 

+3  -  Print  above  plus  the  constraint 
values,  search  direction  vector  and  move  parameter  alpha  on 
each  iteration. 

+4  -  Print  above  plus  gradient  information 
on  each  iteration. 

+5  -  Print  above  plus  each  proposed  design 
vector  and  the  penalized  objective  function  during  the  one- 
dimensional  search  on  each  iteration. 

+  6  -  Debugging  aid  for  optimizer  develop- 
ment.  DESOP  will  pause  after  each  major  operation  is 
performed  during  the  optimization  process. 

DELFUN  -  The  minimum  absolute  change  in  the  objec- 
tive function  to  indicate  convergence  of  the  optimization 
process. 

DABFUN  -  The  minimum  absolute  change  in  the  objec- 
tive function  to  indicate  convergence  of  the  optimization 
process. 

95 


ITMAX  -  The  maximum  number  of  inner  loop  (uncon- 
strained minimizations)  without  increasing  the  penalty 
function. 

ICNDIR  -  The  conjugate  direction  restart  parameter. 
Every  ICNDIR  inner  loop  iterations  a  steepest  descent  search 
is  performed.   It  is  recommended  that  ICNDIR  be  set  equal 
to  NDV  +  1. 

FDCH  -  The  relative  change  in  the  design  variables 
for  calculating  finite  difference  gradients. 

FDCHM  -  The  minimum  absolute  step  in  finite  difference 
gradient  calculations. 

AB0BJ1  -  The  expected  fractional  change  in  the  ob- 
jective function  for  the  first  estimates  of  the  step  size 
to  be  taken  in  the  one-dimensional  search. 

ALPHAX  -  The  maximum  fractional  change  in  any  design 
variable  for  the  first  estimate  of  the  step  size  to  be  taken 
in  the  one-dimensional  search. 

7.   INPUT  CHECK  OF  THE  DESOP  CONTROL  PARAMETERS  -  A  menu 
of  the  optimizer  control  parameters  for  the  DESOP  optimiza- 
tion program  is  displayed.   The  menu  will  display  the  variable 
name,  its  maximum,  minimum,  default  and  present  value.   The 
variables  displayed  are: 

IRMAX  -  The  maximum  number  of  times  that  the  penalty 
parameter  will  be  increased. 

RZ  -  The  starting  value  of  the  penalty  parameter. 
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RMULT  -  The  amount  by  which  ?.Z  is  multiplied  each 

time  ~ha-  it  is  increased. 

EX?G  -   The  amount  by  which  a  vie  lazed  constraint 
value  is  raised  to  an  exponent,  ZX?G. 

0  -  No  scaling  of  z'r.e    desire  variables  is 
perforrr.ee  by  DESO?. 

1  -  The  design  variables  are  scaled  every  ICNIIR 
iteration. 

9.   input  iesii-:;  variables  a:;i  site  CONSTRAINTS  -  A  -e.v: 
of  the  initial  design  variables  and  rjonstrair.es  will  be 
displayed  for  an  existing  analysis  subprogram.   If  this  is 
a  new  analysis  subprcgrarr. ,  the  user  ".•."ill  be  asked  to  enter 
the  initial  values  of   the  design  variables  and  any  side 
constraints  on  the  design  variables. 

A  hard  copy  printout  of  the  optimizer  oontrol  parameters , 
design  variables  and  side  constraints  will  then  be  presented 
to  the  user  to  check  to  ensure  that  the  above  information 
has  been  entered  correctly.   The  user  then  has  the  option 
to  return  to  the  input  routine  and  make  changes  ee  the  above 
information  or  to  continue  and  optimize  his  design.   If  the 
user  chooses  to  continue,  OPCON  will  overlay  the  optimization 
program,  on  CZZZY.    beginning  at  line  C2.   The  analyze  subpro- 
gram, specified  will  then  be  linked  to  the  end  of  the  IESC? 
program.   Ihe  entire  program  will  then  be  STCREI  under  the 
program  name  OP.   Ihe  C?  program,  will  then  be  loaded  with 
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execution  beginning  at  line  0PT1.   The  storing  and  reloading 
process  allows  the  OP  program  to  be  called  and  run  as  a 
separate  program  for  debugging  purposes.   At  the  completion 
of  the  optimization  process,  the  optimizer  program  will 
reload  OPCON  and  return  execution  to  the  OPCON  program.   A 
menu  will  be  displayed  welcoming  the  user  back  to  the  OPCON 
program  giving  him  the  following  options:   to  optimize 
again  using  the  same  ANALIZ  subprogram,  to  optimize  again 
using  a  different  ANALIZ  subprogram,  or  to  terminate  the 
program. 
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APPENDIX  D 
OPCON   Program  Listing 
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APPENDIX  E 
DESOP  Program  Listing 

In  the  development  of  the  DESOP  program,  Refs.  2,  7 ,  8, 
9,  and  10  were  used. 
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APPENDIX  G 
NISCO  Subprogram  Listing 

Below  is  a  cross  reference  list  of  major  equations 
used  in  the  NISCO  program  to  the  references  used  to  develop 
the  equations. 

Reference 

Number       NISCO  Subprogram  Line  Number 

3.  820  -  850 

4.  1085  -  1090 

5.  915,  985  -  1045 

6.  855,  925,  930,  965  -  980,  1100,  1180  - 
1210,  1255  -  1345,  1460,  1465, 

11.  905,  1115,  1125,  1145 

12.  1445  -  1455,  1470  -  1480 

13.  1600  -  1605 
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APPENDIX  H 
Sample  DESOP  Output 

Appendix  H  is  a  sample  computer  output  for  the  DESOP 
program. 
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